Mouse Gastrulation scATAC-seq & scRNA-seq

## Load libraries


suppressPackageStartupMessages(library(ArchR))
library(knitr)
library(rhdf5)
library(uwot)
library(tidyverse)
library(scater)
library(zellkonverter)
library(corrplot)
#library(caret)
h5disableFileLocking()

The scATAC-seq dataset contains two replicates each from 4 stages of mouse gastrulation, E7.5, E8.0, E8.5 and E8.75. We will use the mm10 Genome, because that is the genome.

# names of gastrulation stages
stages <- c("E7.5_rep1", "E7.5_rep2", "E8.0_rep1", "E8.0_rep2", 
 "E8.5_rep1", "E8.5_rep2", "E8.75_rep1", "E8.75_rep2")

# combine all raw fragment files
inputFiles <- setNames(c("raw_files/E7.5_rep1_atac_fragments.tsv.gz",
                       "raw_files/E7.5_rep2_atac_fragments.tsv.gz",
                       "raw_files/E8.0_rep1_atac_fragments.tsv.gz",
                       "raw_files/E8.0_rep2_atac_fragments.tsv.gz",
                       "raw_files/E8.5_rep1_atac_fragments.tsv.gz",
                       "raw_files/E8.5_rep2_atac_fragments.tsv.gz",
                       "raw_files/E8.75_rep1_atac_fragments.tsv.gz",
                       "raw_files/E8.75_rep2_atac_fragments.tsv.gz"),
                       c("E7.5_rep1", "E7.5_rep2", "E8.0_rep1", "E8.0_rep2", 
                         "E8.5_rep1", "E8.5_rep2", "E8.75_rep1", "E8.75_rep2"))

addArchRGenome("mm10")

Create Arrow Files

We will read in the accessible fragments from the fragments file. QC information is automatically computed for each cell. Based on visual inspection of TSS Enrichment score vs log10(nFrags) I used a threshold of at least 1000 unique fragments and a TSS Enrichment score of at least 2.75. At the beginning it is better to be less stringent, because we can still filter the data at a later stage. The cells are filtered based on these threshold and data are stored in arrow files. Genome-wide tile matrices are created using 500-bp bins. The tile matrix T is binary matrix and contains information whether in a particular 500-bp bin/tile i there were any insertions in a particular cell j.

\(T \in \mathbb{R}^{NXD}\) with \(T_{i,j} = 1\) if an insertion is present in tile \(i\) and cell \(j\). \(T_{i,j} = 0\) if no insertion is present in that cell, meaning \(T\) is a binary matrix.

In addition a GeneScoreMatrix is created.

I set the size of the maximum number of fragments included to 1e+07, instead of the default 1e+05, because we cut off some cells otherwise.

ArrowFiles <- createArrowFiles(
  inputFiles = inputFiles,
  sampleNames = names(inputFiles),
  minTSS =2.75, #Dont set this too high because you can always increase later
  minFrags = 1000, # minimum number of mapped fragments required
  # count matrix, instead of using peak it uses fixed-width sliding window of bins across the whole genome
  addTileMat = TRUE, 
  addGeneScoreMat = TRUE, # uses signal proximal to the TSS to estimate gene activity
  subThreading = FALSE,
  maxFrags = 1e+07,
  minFragSize = 10,
  maxFragSize = 2000,
  QCDir = "new_ArchRobject/QualityControl",
  force = FALSE, # change to TRUE if you want to ovewrite the files
  # the length in bp that wraps around nucleosomes -> 
    #identify fragments as sub-nucleosome spanning, mono-nucleosome spanning or multi-nucleosome spanning
  nucLength = 147, 
  # integer vector -> define region up/downstream of TSS to include as promoter region
  # can be used to calculate e.g fraction of reads in promoter (FIP)
  promoterRegion = c(2000, 100),
  # parameters for computing TSS enrichment scores, window (bp) centered at TSS = 101
  # flanking window = 2000 bp
  # norm = size of flank window used for normalization = 100 bp
  # accessibility within 101 bp surrounding the TSS will be normalized to accessibility
  # in 100 bp bins from -2000:-1901bp and 1901: 2000
  TSSParams = list(101, 2000, 100),
  # which chromosomes to exclude form downstream analysis
  # in human and mouse: mitochondrial DNA (chrM) and male sex chromosome (chrY)
  # the fragments are still stored in the arrow files
  excludeChr = c("chrM", "chrY"),
  # number of chunks to divide chromosomes in -> low-memory parallelized reading of input files
  nChunk = 5,
  # name of field in input bam file containing the barcode tag information
  bcTag = "qname",
  offsetPlus = 4, # offset applied to + stranded Tn5 insertion -> account for precise Tn5 binding site
  offsetMinus = -5, 
  logFile = createLogFile("createArrows")
)

Inferring Doublets

Should be removed, because they can interfere with downstream analysis.

Doublet detection-and-removal algorithm: Heterotypic doublets are identified by generating a collection of synthetic doublets. These synthetic doublets are projected into low-dimensional embeddings. Searching for nearest neighbours to the synthetic doublets we can identify doublets in the dataset. This outperforms the prediction of doublets using fragment number (ROC-AUC). (Compared to demuxlet as ground truth)

We can also identify doublets in the scRNA-seq space if we have paired data and remove the cells in this way.

ArrowFiles <- setNames(c("E7.5_rep1.arrow", "E7.5_rep2.arrow", "E8.0_rep1.arrow",
                         "E8.0_rep2.arrow", "E8.5_rep1.arrow", "E8.5_rep2.arrow",
                         "E8.75_rep1.arrow", "E8.75_rep2.arrow"),
                       c("E7.5_rep1", "E7.5_rep2", "E8.0_rep1", "E8.0_rep2", 
                         "E8.5_rep1", "E8.5_rep2", "E8.75_rep1", "E8.75_rep2"))

# for each sample provided doublet information will be assigned to each cell
# this way we can remove doublet-based clusters downstream
doubScores <- addDoubletScores(
  useMatrix = "TileMatrix",
  input = ArrowFiles,
  k = 10, #Refers to how many cells near a "pseudo-doublet" to count.
  nTrial = 5, # number of time to simulate nCell doublets 
  knnMethod = "UMAP", #Refers to the dimensionality reduciton method to use for nearest neighbor search.
  LSIMethod = 1 # oder of normalization: tf-log(idf)
)

Create ArchRProject

An ArchR Project is initialized with some important attributes:

• ouput directory
• sample names
• sampleColData -> matrix containing data for each sample
• cellColData -> contains data associated with each cell
  • after using addDoubletScore() there will be a column for “Doublet Enrichment” and “Doublet Score”
• total number of cells (excluding doublets)
• median TSS score & median number of fragments across all cells and samples

proj <- ArchRProject(
  ArrowFiles = ArrowFiles, 
  outputDirectory = "gastrArchROutput",
  copyArrows = TRUE, #This is recommened so that you maintain an unaltered copy for later usage.
  geneAnnotation = getGeneAnnotation(),
  #genomeAnnotation = getGeneAnnotation(),
  showLogo = FALSE
)

proj

proj <- loadArchRProject("08_ArchRobject_all_timepoints_added_gene_expr/")

# read in existing project, so we don't have to rerun the creation
#saveArchRProject(proj, outputDirectory = "01_FirstQC_ArchRobject/", load = FALSE)

#proj <- loadArchRProject(path = "01_FirstQC_ArchRobject/")

Quality Control

Have a look at this for additional QC! https://bioconductor.org/packages/devel/bioc/vignettes/ATACseqQC/inst/doc/ATACseqQC.html

1. the number of unique nuclear fragments (as opposed to mitochondrial fragments) A cell with very few usable fragments will not provide enough data to make useful conclusions.
2. signal-to-background ratio: if this is low this probably corresponds to dying cells where the entire genome allows random transposition
3. fragment size distribution: since 147 bp are wrapped around a nucleosome it is expected that there are depletions of fragments of this length at regular intervals. We expect to see a periodic distribution of fragmetn size corresponding to nucleosomes (mono, di, tri, …), because Tn5 cannot cut DNA that is tightly wrapped around a nucleosome.

Data before QC and corresponding plots are saved in the Quality Control output folder.

log10(unique fragments) vs TSS enrichment

before QC

We used a TSS enrichment score of > 2.75 and a number of unique fragments > 1000 (log10(1000) = 3). It might make sense to use different thresholds for different samples.

dir <- "QualityControl/"
#all_files <- list.files(path)


beforeQC <- map(stages, function(n) {
  path <- paste0(dir, n, "/")
  all_files <- list.files(path)
  file <- all_files[grepl(paste0(n, "-Pre-Filter-Metadata.rds"), all_files)]
  # read in file for sample n containing original data before QC
  df <- as.data.frame(readRDS(paste0(dir, n, "/", file)))
  df %>% mutate(Sample = n) %>% 
    ggplot(aes(x = log10(nFrags), y = TSSEnrichment)) +
    geom_density2d_filled() +
    geom_hline(yintercept = 2.75, color="green", linetype="dashed") +
    geom_vline(xintercept = log10(1000), color="green", linetype="dashed") +
    #geom_xsidedensity(aes(x=log10(pre_filter_meta$nFrags))) +
    #geom_ysidedensity(aes(y = pre_filter_meta$TSSEnrichment)) +
    ylim(c(2.5, 4.5)) +
    facet_wrap(~Sample) +
    theme(legend.position = "none") +
    labs(x = "Log10 Unique Fragments", y = "TSS Enrichment Score")
})
  
gridExtra::grid.arrange(grobs = beforeQC, ncol = 3)
## Warning: Removed 361 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 382 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 124 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 103 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 75 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 62 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 269 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 277 rows containing non-finite values (stat_density2d_filled).

after QC

metadata <- as.data.frame(getCellColData(proj))
metadata %>% head
##                                 Sample TSSEnrichment ReadsInTSS ReadsInPromoter
## E7.5_rep1#GTATTGCAGCGGATTT-1 E7.5_rep1         3.080       9360           41754
## E7.5_rep1#GCTGTACCACTATGGC-1 E7.5_rep1         3.396      13022           53683
## E7.5_rep1#AGGTTTGGTAGGTTTG-1 E7.5_rep1         3.016       8528           37914
## E7.5_rep1#CACAATATCGGGATTT-1 E7.5_rep1         3.506      11438           45668
## E7.5_rep1#TGTAACTCAAGACTCC-1 E7.5_rep1         3.261      10098           42865
## E7.5_rep1#ATTGGCTAGGCATTAC-1 E7.5_rep1         3.121       7809           34552
##                              ReadsInBlacklist PromoterRatio PassQC
## E7.5_rep1#GTATTGCAGCGGATTT-1             4831     0.1239800      1
## E7.5_rep1#GCTGTACCACTATGGC-1             5458     0.1661714      1
## E7.5_rep1#AGGTTTGGTAGGTTTG-1             5621     0.1282907      1
## E7.5_rep1#CACAATATCGGGATTT-1             5072     0.1627315      1
## E7.5_rep1#TGTAACTCAAGACTCC-1             4411     0.1639598      1
## E7.5_rep1#ATTGGCTAGGCATTAC-1             5084     0.1334044      1
##                              NucleosomeRatio nMultiFrags nMonoFrags nFrags
## E7.5_rep1#GTATTGCAGCGGATTT-1        1.368854       34465      71085 168390
## E7.5_rep1#GCTGTACCACTATGGC-1        1.327239       32948      69408 161529
## E7.5_rep1#AGGTTTGGTAGGTTTG-1        1.485760       32186      59445 147766
## E7.5_rep1#CACAATATCGGGATTT-1        1.396860       28613      58542 140317
## E7.5_rep1#TGTAACTCAAGACTCC-1        1.316832       27106      56421 130718
## E7.5_rep1#ATTGGCTAGGCATTAC-1        2.599450       28797      35978 129501
##                              nDiFrags DoubletScore DoubletEnrichment
## E7.5_rep1#GTATTGCAGCGGATTT-1    62840    138.38655             6.700
## E7.5_rep1#GCTGTACCACTATGGC-1    59173      1.53785             1.750
## E7.5_rep1#AGGTTTGGTAGGTTTG-1    56135     14.75493             2.600
## E7.5_rep1#CACAATATCGGGATTT-1    53162    323.30622            10.925
## E7.5_rep1#TGTAACTCAAGACTCC-1    47191      0.00000             0.400
## E7.5_rep1#ATTGGCTAGGCATTAC-1    64726     67.55402             4.650
##                              BlacklistRatio            celltypes Gex_nUMI
## E7.5_rep1#GTATTGCAGCGGATTT-1     0.01434468                  Gut     5872
## E7.5_rep1#GCTGTACCACTATGGC-1     0.01689480                  Gut    14826
## E7.5_rep1#AGGTTTGGTAGGTTTG-1     0.01901994 Rostral_neurectoderm     3992
## E7.5_rep1#CACAATATCGGGATTT-1     0.01807336         ExE_ectoderm    19045
## E7.5_rep1#TGTAACTCAAGACTCC-1     0.01687220             Epiblast    12849
## E7.5_rep1#ATTGGCTAGGCATTAC-1     0.01962919         ExE_ectoderm     3824
##                              Gex_nGenes Gex_MitoRatio Gex_RiboRatio
## E7.5_rep1#GTATTGCAGCGGATTT-1       3328             0             0
## E7.5_rep1#GCTGTACCACTATGGC-1       5295             0             0
## E7.5_rep1#AGGTTTGGTAGGTTTG-1       2508             0             0
## E7.5_rep1#CACAATATCGGGATTT-1       6189             0             0
## E7.5_rep1#TGTAACTCAAGACTCC-1       4926             0             0
## E7.5_rep1#ATTGGCTAGGCATTAC-1       2319             0             0

metadata %>% group_by(Sample) %>% summarise(n = n())
## # A tibble: 8 × 2
##   Sample         n
##   <chr>      <int>
## 1 E7.5_rep1   6675
## 2 E7.5_rep2   7955
## 3 E8.0_rep1   5154
## 4 E8.0_rep2   4526
## 5 E8.5_rep1   7702
## 6 E8.5_rep2   6310
## 7 E8.75_rep1  3279
## 8 E8.75_rep2  4385

metadata %>% ggplot() +
  geom_density2d_filled(aes(x = log10(nFrags), y = TSSEnrichment)) +
  geom_hline(yintercept = 2.75, color="green", linetype="dashed") +
  geom_vline(xintercept = log10(1000), color="green", linetype="dashed") +
  #geom_xsidedensity(aes(x=log10(pre_filter_meta$nFrags))) +
  #geom_ysidedensity(aes(y = pre_filter_meta$TSSEnrichment)) +
  ylim(c(2.5, 4.5)) +
  facet_wrap(~Sample) +
  theme(legend.position = "none") +
  labs(x = "Log10 Unique Fragments", y = "TSS Enrichment Score")
## Warning: Removed 132 rows containing non-finite values (stat_density2d_filled).

ArchR plots

```#{r, fig.width=8, fig.height=15} #create 3 separate dataframes for all samples samples <- map(unique(proj\(Sample), function(name){ index <- BiocGenerics::which(proj\)Sample %in% name) cells <- proj$cellNames[index] sample_subset <- proj[cells] df <- getCellColData(sample_subset, select = c(“log10(nFrags)”, “TSSEnrichment”)) p <- ggPoint( x = df[, 1], y = df[, 2], colorDensity = TRUE, # should the density of points on the plot be indicated by color? continuousSet = “sambaNight”, xlabel = “Log10 unique fragments”, ylabel = “TSS enrichment”, title = paste0(“Sample:”, name), #xlim = c(log10(500), quantile(df[,1], probs = 0.99)), #ylim = c(0, quantile(df[,2], probs = 0.99)) ) + geom_hline(yintercept = 2.75, lty = “dashed”) + geom_vline(xintercept = 3, lty = “dashed”) })

#do.call(“grid.arrange”, c(samples, ncol=4))

gridExtra::grid.arrange(grobs = samples, ncol = 2)

### Plotting sample statistics

* when we have distinct samples, it can be important to compare various
metric between samples
* ridge plots & violin plots are used for grouped data


```r

p1 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_")) %>% 
  ggplot() +
  geom_density(aes(x = TSSEnrichment, fill = Sample), alpha = 0.8) 
## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

p2 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_")) %>% 
  ggplot() +
  ggridges::geom_density_ridges(aes(x = TSSEnrichment, y = Sample,
                                    fill = Sample), alpha = 0.8) +
  theme(legend.position = "none")

p3 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_")) %>% 
  ggplot() +
  geom_violin(aes(x = Sample, y = TSSEnrichment, fill = Sample), alpha = 0.8) +
  geom_boxplot(aes(x = Sample, y = TSSEnrichment,fill = Sample), alpha = 0.4) + 
  theme(legend.position = "none") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none") +
  labs(title = "TSS Enrichment")


cowplot::plot_grid(p3, p2, p1, ncol = 3)

p1 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_"), log10_nFrags = log10(nFrags)) %>% 
  ggplot() +
  geom_density(aes(x = log10_nFrags, fill = Sample), alpha = 0.8)

p2 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_"), log10_nFrags = log10(nFrags)) %>% 
  ggplot() +
  ggridges::geom_density_ridges(aes(x = log10_nFrags, y = Sample,
                                    fill = Sample), alpha = 0.8) +
  theme(legend.position = "none") 


p3 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_"), log10_nFrags = log10(nFrags)) %>% 
  ggplot() +
  geom_violin(aes(x = Sample, y = log10_nFrags, fill = Sample), alpha = 0.8) +
  geom_boxplot(aes(x = Sample, y = log10_nFrags,fill = Sample), alpha = 0.4) + 
  theme(legend.position = "none") +
  labs(title = "number of fragments") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none") 

cowplot::plot_grid(p3, p2, p1, ncol = 3)

Plot Fragment Size Distribution & TSS Enrichment Profiles

• the distribution of fragments size can be very different between samples, cell types and batches -> these differences do not necessarily correlate with differences in quality
• the dip is the fragment size of a nucleosome ~147bp
• TSS enrichment profiles
  • clear peak in the center
  • smaller shoulder peak right of the center caused by well positioned +1 nucleosome

#{r, fig.width=5, fig.height=5} p1 <- plotFragmentSizes(ArchRProj = proj) p2 <- plotTSSEnrichment(ArchRProj = proj) ggAlignPlots(p1, p2, type = "v")

Filtering Doublets

In the plots below the Doublet Enrichment represents the enrichment of simulated dobulets nearby each cell, compared to the expected number of doublets if a uniform distriution is assumed.

The doublet density vissualizes where the synthetic doublets are located after projecting them into the 2D embedding.

TODO: Chage path and directory to different name, because the are functions!

# the doublet information is saved in a simpleListobject



dir <- "QualityControl/"


doublet_enrichment <- map(stages, function(n) {
  path <- paste0(dir, n, "/")
  all_files <- list.files(path)
  file <- all_files[grepl(paste0(n, "-Doublet-Summary.rds"), all_files)]
  # read in file containing doublet summary for sample n 
  doublet_summary <- readRDS(paste0(dir, n, "/", file))
  doublet_summary[[1]] %>% 
  ggplot() +
  geom_point(aes(x = X1, y = X2, col = enrichment), size = .8) +
  scale_color_viridis_c() +
  guides(fill=guide_legend(title="DoubletEnrichment")) +
  labs(title = paste0("Simulated Doublet Enrichment", n))
})
  
gridExtra::grid.arrange(grobs = doublet_enrichment, ncol = 2)

doublet_density <- map(stages, function(n) {
  path <- paste0(dir, n, "/")
  all_files <- list.files(path)
  file <- all_files[grepl(paste0(n, "-Doublet-Summary.rds"), all_files)]
  # read in file containing doublet summary for sample n 
  doublet_summary <- readRDS(paste0(dir, n, "/", file))
  doublet_summary[[2]] %>% 
  ggplot() +
  geom_point(aes(x = x, y = y, col = density), size = .1) +
  geom_point(data = doublet_summary[[1]], aes(x = X1, y = X2), size = .1,
             alpha = .3) +
  scale_color_viridis_c() +
  guides(fill=guide_legend(title="Simulated Doublet Density")) +
  labs(title = "Simulated doublet density overlayed")

})

gridExtra::grid.arrange(grobs = doublet_density, ncol = 2)

With the function addDoubleScores() information on predicted doublets has been added. Now we will filter the putative doublets. ArchR automatically prints the number of cells removed from each sample and the corresponding percentage.

arguments:

• cutEnrich = minimum cutoff for DoubletEnrichment, number of simulated doublets divided by expected number given a random uniform distribution
• cutScore = minimum cutoff for Doublet Score, represents -log10(binomial adjusted p-value) for the DoubletEnrichmentadd
• filterRatio = maximum ratio of predicted doublets to filter based on number of pass-filter cells (A higher filterRatio means that more cells are removed) e.g. 5000 cells

maximum would be filterRatio * 5000 / 100000 = filterRatio * 5000 * 0.05

This way samples with different percentage of doublets will be filtered accordingly.

As can be seen below, in E7.5_rep1 no doublets are filtered. Why might that be?

# in our case we now have 10 251 cells as opposed to 10 661 cells before
# filtering -> 410 cells were removed (3.85%)
proj <- filterDoublets(ArchRProj = proj)

Use annotated cells

Subset cells

What is the overlap between my QC and the processed and annotated data?

# read in processed scATAC-seq data
proc_atac <- readH5AD("jupyter_notebooks/anndata_atac_peak_matrix.h5ad")
proc_rna <- readH5AD("jupyter_notebooks/anndata_rna.h5ad")

print(paste0("The number of genes in the processed data is: ", dim(proc_rna)[[1]]))
## [1] "The number of genes in the processed data is: 32285"
print(paste0("The number of peaks in the processed data is: ", dim(proc_atac)[[1]]))
## [1] "The number of peaks in the processed data is: 180499"
print(paste0("The number of cells in the processed data is: ", dim(proc_atac)[[2]]))
## [1] "The number of cells in the processed data is: 45991"
print(paste0("The number of cells in the ArchR object is: ", dim(proj)[[2]]))
## [1] "The number of cells in the ArchR object is: "

overlap <- proj$cellNames[proj$cellNames %in% colnames(proc_atac)]
print(paste0("Out of the ", dim(proc_atac) [[2]], " cells found in the processed data, ", length(overlap), " are also found in our ArchR object"))
## [1] "Out of the 45991 cells found in the processed data, 45986 are also found in our ArchR object"

filt_proj <- proj[proj$cellNames[proj$cellNames %in% colnames(proc_atac)], ]

print(paste0("The number of remaining cells is: ", dim(filt_proj)[[2]]))
## [1] "The number of remaining cells is: "

filt_proc_rna <- proc_rna[ , filt_proj$cellNames]
filt_proc_atac <- proc_atac[ , filt_proj$cellNames]

gene_scores <- getMatrixFromProject(filt_proj, useMatrix = "GeneScoreMatrix")

#filt_proj <- loadArchRProject(path = "05_newdata_ArchRObject/")

Save the subset ArchR project and proceed with it for now. It might be interesting to revisit and see if I can cluster and identify the cell types de novo at a later stage!

Dimensionality reduction

Because we can have maximally two accessible alleles per cell, the scATAC-seq data is sparse. Therefore, the majority of accessible regions are not transposed, meaning that most loci will have 0 accessible alleles. A zero could mean “non-accessible” or “not sampled”. For many analysis we can use a binarized matix. Imporantly, the 1s have information, BUT the 0s do not!

A PCA would result in high inter-cell similarity at all 0 positions. An alternative approach for dimensionality reduction is a layered dimensionality reduction. First, Latent Semantic Indexing (LSI) is used. LSI is an approach from language processing. Different samples are the “documents” and different regions/peaks are the “words”.

Iterative LSI

1. compute term frequency (depth normalization to 10,000 per single cell)

\(TF = \frac{C_{ij}}{F_{j}}\), the term frequency of \(C_{ij}\) which is tiles i in cell j and \(F_{j}\) being the total number of tiles with insertions in cell j.

2. Inverse document frequency

• weights features by how often they occur
• more weight to less frequent insertions

\(IDF_{i} =\log {\frac{N}{1 + DF_{i}}}\) with N being the total number of cells in the dataset and \(DF_{i}\) being the total number of cells (documents) in which the tile i is present.

3. The term frequency TF is normalized by the inverse document frequncy IDF. You get a \(TF-IDF\) matrix (term frequency-inverse document frequency) which tells us how important a region/peak is to a sample. In other words you transform a binary matrix to a non-binary matrix.

\(TF-IDF = \log{1 + TF{ij}} * IDF_{i} * 10^{4}\)

4. SVD identifies the most valuable information across samples. Then we can use these most valuable features to represent the data in a lower dimensional space
5. Clusters are identified with Seurat’s Shared Nearest Neighbor clustering
6. Sum accessibility across all single cells in each cluster -> log-normalize
7. Identify most variable features across the clusters
8. repeat with most variable peaks as features

With LSI we can reduce the dimensionality of the sparse insertion matrix to tens or hundreds. Then UMAP or t-SNE can be used to visualize the data

Unlike in scRNA-seq we cannot select the top highly variable features before dimensionality reduction (high noise, low reproducibility). Rather the iterative LSI approach first computes a LSI on the most accessible tiles (this will identify clusters corresponding to the major cell types). Then, ArchR computes the average accessibility across these clusters across all features. Next, the most variable peaks across these clusters are identified. The most highly accessible peaks are the features of a new round of LSI. We can set how many rounds of LSI we want to be peformed.

Using iterative LSI reduces batch effects. If you see some batch effects you could try to add more LSI iterations and start from a lower initial clustering resolution. Also, the number of variable features can be lowered. #

#filt_proj <- addIterativeLSI(ArchRProj = filt_proj, useMatrix = "TileMatrix", 
                             #name = "atac_LSI_25000")
#filt_proj <- addIterativeLSI(ArchRProj = filt_proj, useMatrix = "TileMatrix", 
                             #name = "atac_LSI_50000")
filt_proj <- addIterativeLSI(ArchRProj = filt_proj, useMatrix = "TileMatrix", 
                            name = "atac_LSI_100000")

Clustering

Here we do not need to perform clustering, because we already have celltype annotations.

Visualization

We will add the celltype annotation from the processed data to our ArchR project. Since we already have annotations, we do not need to perform clustering.

#filt_proj <- addUMAP(ArchRProj = filt_proj, reducedDims = "atac_LSI_25000", force = TRUE)

#filt_proj <- addUMAP(ArchRProj = filt_proj, reducedDims = "atac_LSI_50000", force = TRUE)

filt_proj <- addUMAP(ArchRProj = filt_proj, reducedDims = "atac_LSI_100000", 
                     force = TRUE, verbose = FALSE)

Add the celltype annotations to the ArchRobject

df <- colData(filt_proc_atac) %>% as.data.frame() %>% rownames_to_column("cell") #%>%
  #filter(cell %in% filt_proj$cellNames)

filt_proj <- addCellColData(
  ArchRProj = filt_proj,
  data = setNames(as.character(df$celltype.mapped_seurat), as.character(df$cell)),
  name = "celltypes",
  cells = df$cell
)

# add custom colors


colPalette_celltypes = c('#532C8A',
 '#c19f70',
 '#f9decf',
 '#c9a997',
 '#B51D8D',
 '#3F84AA',
 '#9e6762',
 '#354E23',
 '#F397C0',
 '#ff891c',
 '#635547',
 '#C72228',
 '#f79083',
 '#EF4E22',
 '#989898',
 '#7F6874',
 '#8870ad',
 '#647a4f',
 '#EF5A9D',
 '#FBBE92',
 '#139992',
 '#cc7818',
 '#DFCDE4',
 '#8EC792',
 '#C594BF',
 '#C3C388',
 '#0F4A9C',
 '#FACB12',
 '#8DB5CE',
 '#1A1A1A',
 '#C9EBFB',
 '#DABE99',
 '#65A83E',
 '#005579',
 '#CDE088',
 '#f7f79e',
 '#F6BFCB')

celltypes <- (as.data.frame(getCellColData(filt_proj)) %>% group_by(celltypes) %>% 
 summarise(n = n()))$celltypes

col <- setNames(colPalette_celltypes, celltypes)

df <- as_data_frame(cbind(getCellColData(filt_proj), getEmbedding(filt_proj)) ) %>%
  rename(c(umap1 = atac_LSI_100000.UMAP_Dimension_1, umap2  = atac_LSI_100000.UMAP_Dimension_2))

variables <- c("celltypes", "Sample", "nFrags", "DoubletScore")


plot1 <- map(c("celltypes"), function(n){
  ggplot() +
  geom_point(aes(x = df %>% pull("umap1"), y = df %>% pull("umap2"), 
                 col = df %>% pull(n)),  size = .8) +
    guides(col=guide_legend(title=paste0(n))) +
    scale_color_manual(values = col)+
    xlab("umap1") +
    ylab("umpa2") +
    labs(title = paste0(n))
})

plots2 <- map(c("nFrags", "DoubletScore"), function(n){
  ggplot() +
  geom_point(aes(x = df %>% pull("umap1"), y = df %>% pull("umap2"), 
                 col = df %>% pull(n)), size = .8) +
    scale_color_viridis_c() +
    guides(fill=guide_legend(title=paste0(n))) +
    xlab("umap1") +
    ylab("umpa2") +
    labs(title = paste0(n))
})

plot3 <- map(c("Sample"), function(n){
  ggplot() +
  geom_point(aes(x = df %>% pull("umap1"), y = df %>% pull("umap2"), 
                 col = df %>% pull(n)),  size = .8) +
    guides(col=guide_legend(title=paste0(n))) +
    xlab("umap1") +
    ylab("umpa2") +
    labs(title = paste0(n))
})


plot1
## [[1]]

do.call(gridExtra::grid.arrange, c(plot3, plots2, ncol=3))#, nrow = 2))

#saveArchRProject(filt_proj, outputDirectory = "02_filtered_annotated_ArchRobject", load = FALSE)
#saveArchRProject(filt_proj, outputDirectory = "05_newdata_ArchRObject/", load = FALSE)

TODO: Continue here!

Check gene scores for the different celltypes

filt_proj <- addImputeWeights(filt_proj, reducedDims = "atac_LSI_100000")
## Warning in sprintf("Completed Getting Magic Weights!",
## round(object.size(weightList)/10^9, : one argument not used by format 'Completed
## Getting Magic Weights!'

markerGenes  <- c(
    "Dnmt3b",  # Epigblas
    "Ptn", # Forebrain/Midbrain/Hindbrain
    "Pou5f1", # Primitive Streak
    "Myl7", # cardiomyocytes
    "Hoxa10" #Allantois

  )

#{r, fig.height=10, fig.width=10} p <- plotEmbedding( ArchRProj = filt_proj, colorBy = "GeneScoreMatrix", name = markerGenes, embedding = "UMAP", imputeWeights = getImputeWeights(proj) ) do.call(gridExtra::grid.arrange, c(p, ncol = 3))

Visualizing Genome Browser Tracks

Browse local chromatin accessibility at marker genes. Plot genome browser tracks per celltype

p <- plotBrowserTrack(
    ArchRProj = filt_proj, 
    groupBy = "celltypes", 
    geneSymbol = markerGenes, # the plot window is centered at the TSS
    upstream = 50000,
   
     downstream = 50000
)

• Dnmt3b = Epiblast
• Ptn = Forebrain/Midbrain/Hindbrain
• Pouf5f1 = Primitive Streak
• Myl7 = cardiomyocytes
• Hoxa10 = Allantois

Since we group the cells accoring to celltype below and plot the average signal/accessibility. Since some signals are only expected in a subset of the embryonic stages, the signals might be diluted if we look at all stages at once and take the mean accessibility over different stages. It would probabl

```#{r,fig.height=15, fig.width=10} for (elem in p){ grid::grid.newpage() grid::grid.draw(elem) }

# Different Time points

To differentiate between expression at different stages we will have to split
the ArchR projects into 

TODO

# Combine with RNA-seq data

Since we are using multiome data, we do not need to integrate the ATAC and RNA
data, but we have both informations from the same cell. We will use the processed
scRNA-seq data which we filtered to contain only the cells which overlap with 
the cells in the ArchR object.

Look here for more details on multiome analysis in ArchR:
https://greenleaflab.github.io/ArchR_2020/Ex-Analyze-Multiome.html

One interesting thing to notice is that the number of genes in the gene score 
matrix is different from the number of genes in the RNA dataset.

I get an error when trying to extract the tile matrix. Not sure what is happening
there. 


```r
gene_score_mat <- assays(gene_scores)[[1]]

# set names of the gene score matrix
rownames(gene_score_mat) <- rowData(gene_scores)$name

print(paste0("The number of genes in the ATAC-seq Gene Score Matrix is: ",
             nrow(gene_score_mat)))
## [1] "The number of genes in the ATAC-seq Gene Score Matrix is: 24333"

print(paste0("The number of genes in the RNA-seq Gene Expression Matrix is: ",
             dim(filt_proc_rna)[1]))
## [1] "The number of genes in the RNA-seq Gene Expression Matrix is: 32285"

print(paste0("The number of genes found in both the Gene Expression Matrix and the Gene Score Matrix:", 
      length(intersect(rownames(gene_score_mat), rownames(filt_proc_rna)))))
## [1] "The number of genes found in both the Gene Expression Matrix and the Gene Score Matrix:22273"

Filter genes which are found in gene score matrix and gene expression matrix. Reorder cell names to match between gene scores and gene expression. Here it would be nice to show VennDiagrams to quickly be able to understand what the differences in genes are between ArchR object and gene expression matrix. Is there a way to include also genes which do not overlap?

# Try to reorder the cellnames to the 

# this is the order of the cell names in the ArchR project gene score matrix
order <- rownames(getCellColData(filt_proj))

# extract the gene expression matrix 
rna_counts <- assays(filt_proc_rna)[[1]]

# reorder gene expression matrix cell names to match the cellnames of the ArchR object
rna_counts <- rna_counts[, order]

# I checked that they are the same
#colnames(rna_counts) == rownames(getCellColData(filt_proj))

# get the gene names of the gene score matrix
genes <- rownames(gene_score_mat)

# subset the gene expression matrix to obtain only the genes found in the gene 
# score matrix
rna_counts <- rna_counts[rownames(rna_counts) %in% genes,]
gene_score_mat <- gene_score_mat[rownames(gene_score_mat) %in% rownames(rna_counts),]


rm(genes)
rm(proc_atac)
rm(proc_rna)

Next we have to create a GRanges object which will be used in the next step to create a RangedSummarizedExperiment. When using the gene annotation metadata stored in the rowData of the RNA-SCE for creating this RangedSummarizedExperiment with the ArchRobject, I got an error message. This might have been due to the fact that the gene annotations looked different between the RNA-SCE and the ArchR object. Therefore, I tried using the same gene annotation as ArchR uses. In order for this to work, we will have to exclude any genes which are not used by ArchR.

# get the Granges object from the Archr project
archr_granges <-getGenes(filt_proj)

# filter for genes which are included in the gene expression matrix
granges <- subset(archr_granges, symbol %in% rownames(rna_counts))

# the problem here was that the gene annotations of the expression matrix
did not match the gene annotations of the Archr object.

df <- rowData(rna_overlap) %>% as.data.frame() %>% filter(gene %in%genes)
df

granges <- df %>% GRanges()
granges

We will now add the gene expression matrix to the ArchR object. First, we have to convert the Single Cell Experiment to a Summarized Experiment.

# create summarized experiment of the scRNA-seq
seRNA <- SummarizedExperiment(assays = list(counts=rna_counts),
                          colData = colnames(filt_proc_rna),
                          rowRanges = granges)
seRNA


# add RNA gene expression matrix to the ArchR object
#filt_proj <- addGeneExpressionMatrix(input = filt_proj, seRNA = seRNA, force = TRUE)


# save the ArchR object with the added gene expression matrix
#saveArchRProject(filt_proj, outputDirectory = "08_ArchRobject_all_timepoints_added_gene_expr/", load = FALSE)

We will add a second LSI for the RNA data. The one based on the ATAC data is called “LSI_ATAC”. We will call the one for the RNA data “LSI_RNA”

marker_genes <- c("Lamb1", "Plat", "Sparc", "Elf5", "Ascl2", "Tfap2c", "Ttr",
                  "Apoa2", "Apoe", "Cystm1", "Emb", "Spink1", "Amot", "Krt19")

p <- plotEmbedding(
    ArchRProj = filt_proj, 
    colorBy = "GeneScoreMatrix", 
    name = marker_genes, 
    embedding = "UMAP",
    quantCut = c(0.01, 0.95),
    imputeWeights = NULL
)

p <- plotEmbedding(
    ArchRProj = filt_proj, 
    colorBy = "GeneScoreMatrix", 
    name = marker_genes, 
    embedding = "UMAP",
    quantCut = c(0.01, 0.95),
    imputeWeights = NULL
)

Visualize the gene activity scores

heatmap(as.matrix(gene_score_mat[1:5, 1:5]))

print(paste0("Minimum gene score: ", min(gene_score_mat)))
## [1] "Minimum gene score: 0"
print(paste0("Maximum gene score: ", max(gene_score_mat)))
## [1] "Maximum gene score: 237.301"

print(paste0("Number of genes in the gene activity score matrix: ", dim(gene_score_mat)[[1]]))
## [1] "Number of genes in the gene activity score matrix: 22273"

tibble(gene = rownames(gene_score_mat),
          mean = Matrix::rowMeans(gene_score_mat),
          var = MatrixGenerics::rowVars(gene_score_mat)) %>% 
  mutate(log10_mean = log10(mean),
         log10_var = log10(var)) %>% ggplot() +
  geom_point(aes(x = log10_mean, y = log10_var)) +
  labs(title = "Gene activity score from scATAC-seq")

Visualize the raw gene expression counts

print(paste0("Number of genes in gene expression matrix: ", dim(rna_counts)[[1]]))
## [1] "Number of genes in gene expression matrix: 22273"

tibble(gene = rownames(rna_counts),
          mean = Matrix::rowMeans(rna_counts),
          var = MatrixGenerics::rowVars(rna_counts)) %>% 
  mutate(log10_mean = log10(mean),
         log10_var = log10(var)) %>% ggplot() +
  geom_point(aes(x = log10_mean, y = log10_var)) +
  labs(title = "Gene expression raw counts")

Normalize the RNA counts using log2 transformation:

1. normalize by library size
2. scale with factor 1000
3. add pseudocount
4. log transform

norm_rna <- Matrix::t(Matrix::t(rna_counts) /
                        Matrix::colSums(rna_counts))

#scaling factor = 1000
sc_factor <- 1e4

# add pseudocount of 1
lognorm_rna <- (log(norm_rna * sc_factor + 1))

tibble(gene = rownames(lognorm_rna),
          mean = Matrix::rowMeans(lognorm_rna),
          var = rowVars(lognorm_rna)) %>% 
  mutate(log10_mean = log10(mean),
         log10_var = log10(var)) %>% ggplot() +
  geom_point(aes(x = log10_mean, y = log10_var)) +
  labs(title = "Gene expression log-normalized RNA expression")

Correlation between gene activity score and gene expression:

marker_genes <- c("Lamb1", "Plat", "Sparc", "Elf5", "Ascl2", "Tfap2c", "Ttr",
                  "Apoa2", "Apoe", "Cystm1", "Emb", "Spink1", "Amot", "Krt19")

```#{r, fig.width=10, fig.height=15} plots <- map(marker_genes, function(n){ plot <- tibble(gene_expr = lognorm_rna[n, ], gene_score = gene_score_mat[n, ]) %>% ggplot(aes( x = gene_expr,y = gene_score)) + geom_point() + geom_smooth(formula = y ~ x, method = “lm”, size=0.1)

})

gridExtra::grid.arrange(grobs = plots, ncol = 3)

test <- cor(rna_counts[“Apoe”, ], gene_score_mat[“Lamb1”, ]) test

```#{r}
plot <- tibble(gene_expr = lognorm_rna["Lamb1", ], 
       gene_score = gene_score_mat["Lamb1", ]) %>% 
  ggplot(aes( x = gene_expr,y = gene_score)) + geom_point() +
  geom_smooth(formula = y ~ x, method = "lm", size=0.1) 
  #ggpubr::stat_cor(aes(label = ..r.label..), method = "pearson", r.digits=2)
  
plot

test <- ArchR::getMatrixFromProject(filt_proj, useMatrix = "GeneExpressionMatrix")

colData(test) %>% as.data.frame %>% colnames() %>% tail(1) %>% print()

filt_proj <- addIterativeLSI(
    ArchRProj = filt_proj, 
    clusterParams = list(
      resolution = 0.2, 
      sampleCells = 10000,
      n.start = 10
    ),
    saveIterations = FALSE,
    useMatrix = "GeneExpressionMatrix", 
    depthCol = "Gex_nUMI",
    varFeatures = 2500,
    firstSelection = "variable",
    binarize = FALSE,
    name = "LSI_RNA"
)
ArchR::getCellColData(filt_proj) %>% as.data.frame()

filt_proj <- addIterativeLSI(
  ArchRProj = filt_proj, 
  saveIterations = FALSE,
  clusterParams  = list(
    resolution = 0.2,
    sampleCells = 10000,
    n.start = 10
  ),
  useMatrix = "GeneExpressionMatrix",
  depthCol = "Gex_nUMI",
  varFeatures = 100000, 
  firstSelection = "Var", 
  binarize = FALSE,
  name = "LSI_RNA"
)

Now we will combine the two embeddings.

#Combined Dims
proj <- addCombinedDims(proj, reducedDims = c("LSI_ATAC", "LSI_RNA"), name =  "LSI_Combined")

#UMAPs
proj <- addUMAP(proj, reducedDims = "LSI_ATAC", name = "UMAP_ATAC", minDist = 0.8, force = TRUE)

proj <- addUMAP(proj, reducedDims = "LSI_RNA", name = "UMAP_RNA", minDist = 0.8, force = TRUE)

proj <- addUMAP(proj, reducedDims = "LSI_Combined", name = "UMAP_Combined", minDist = 0.8, force = TRUE)